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ABSTRACT 

We introduce a new maximum likelihood method to model the density profile of Blue Hor- 
izontal Branch and Blue Straggler stars and apply it to the Sloan Digital Sky Survey Data 
Release 8 (DR8) photometric catalogue. There are a large number 20,000) of these trac- 
ers available over an impressive 14, 000 deg^ in both Northern and Southern Galactic hemi- 
spheres, and they provide a robust measurement of the shape of the Milky Way stellar halo. 
After masking out stars in the vicinity of the Virgo Overdensity and the Sagittarius stream, 
the data are consistent with a smooth, oblate stellar halo with a density that follows a broken 
power-law. The best fitting model has an inner slope a-^^i ~ 2.3 and an outer slope ctout ^ 4.6, 
together with a break radius occurring at ^ 27 kpc and a constant halo flattening (that is, ratio 
of minor axis to major axis) of q ^ 0.6. Although a broken power-law describes the density 
fall-off most adequately, it is also well fit by an Einasto profile. There is no strong evidence 
for variations in flattening with radius, or for triaxiality of the stellar halo. 

Key words: galaxies: general - galaxies: haloes - galaxies: individual: Milky Way - Galaxy: 
stellar content - Galaxy: structure - galaxies: photometry 



1 INTRODUCTION 

The time required for stars in the stellar halo to exchange their en- 
ergy and angular momentum is very long compared to the age of 
the Galaxy. Therefore, such stars preserve memories of their ini- 
tial conditions, and so the structure of the stellar halo is intimately 
linked to the formation mechanism of the Galaxy i tself. This fun- 
damental insight was already noted bv lEggen et al.l ( Il962[) . It is the 
reason why the stellar halo has attracted such interest despite con- 
taining only a small fraction of the total stellar mass of the Galaxy. 
■ Stars diffuse more quickly in configuration space, as opposed to en- 
ergy and angular momentum space. So, the spatial structure of the 
stellar halo may be smooth, even though it is built up from merging 
and accretion. 

The simplest way of studying the stellar halo is through star- 
counts. Typically, RR Lyrae or blue horizontal branch stars (BHBs) 
are used as tracers, a s they are relatively bright {Mg 0.5 — 0.7, 
e.g. lSirko et al .120041) and can be detected at radii out to ~ 100 kpc. 
The gathering of such data is painstaking work, and carries the price 
that sample sizes are often small. Such studies are consistent with a 
stellar halo that is round in the outskirts (with a minor-to- major axis 
ratio q = 1) and more flattened in th e inner parts with g 0.5 (e.g., 
lHartwick|[l987l: l Preston et al. Rather than selecting typical 

halo stars, an alternative approach is to model deep star count data 
in pencil-beam surveys at intermediate and high galactic latitudes, 
allowing for contamination of the starcounts by the thin and thick 
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disk populations. This was attempted bv lRobin et al] bOOOh . who 
found a best-fit halo density law with flattening q ^ 0.76, together 
with a power-law fall-off a of 2.4 (t hat is, p ^ (dis tance)""). A 
similar, slightly later, attempt by iSiegel et al. I( l2002l) using data in 
seven Kapteyn selected areas yielded q ~ 0.6 and a ~ 2.75. 

Efforts to detect variations in the flat tening with radius have 
also been undertaken. IPreston et al.l ( Il99lh argued that the flatten- 
ing changes from strongly flattened {q — 0.5) at 1 kpc to almost 
round at 20 kpc. However, work by Isiuis & Arnold! ( Il998l) using 
a compilation of 340 RR Lyraes and BHBs found a constant flat- 
tening of g ~ 0.5 with no evidence for changes with radius, as 
we ll as and a power-law in dex a ~ —3.2. The most recent work 
by iDe Propris et al.l JloTol) utilising 666 BHB stars from the 2dF 
Quasar Redshift Survey find that the halo is approximately spheri- 
cal with a power-l aw index of « = 2.5 out to ~ 100 kpc. Similarly, 
ISesaret al.l ( l201lh studying Main Sequence Tum-Off (MSTO) stars 
from the Canada-France-Hawaii Telescope Legacy Survey find that 
the flattening is approximately constant at g ~ 0.7 out to 35 kpc. 

The Sloan Digital Sky Survey (SDSS) has transformed our 
knowledge of the stellar halo. Although it had been suspected 
that the stellar halo is criss-crossed with streams and substruc- 
tures ever since the discovery of the disrupting Sagittarius (Sgr), 
the SDSS provid ed a memorable pictu re of the debris in The 
Field of Streams (iBelokurov et al. I l2006h . A wealth of substruc- 
ture has now been identified, including the Sagittarius stream, the 



Virgo Overdensity and t he He r cules-Aquila Cloud (e.g. Ilbata et al 



19951: IBelokurov et"ai] l2006l : ljuric et al.l l2008l : IBelokurov et~al 



2007). This has been seen as vindication of modem theories of 
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galaxy formation, which predict that stellar haloes are built up 
almost exclusively from t he debris of disrupting s atellites (e.g. 
Bullock & Johnstoij|200i lOe Lucia & HelmJl2008l : ICooper et al.] 



201Cl) . A number of studies have attempted to model the smooth 



halo component by avoiding these known substructures (e.g. 
ljuric et ^ [ 2008h . The results are only in very rough agreement, 
suggesting that the density profile of the Milky Way has a power- 
law slope in the range 2 < a < 4 and a flattening varying 
from 0.4 < q < 0.8 ( e.g. |Y annv et al.' 2000; Ch en et al.|[200lf ; 
iNewberg & Yannvll2006l ; ljurig et al..2008; Sesar et al.bOl ih . 

Even though panoramic photometric surveys like SDSS do 
provide a large sample of stellar halo tracers over a considerable 
portion of the sky, there is no consensus on the flattening and shape 
of the stellar halo. MSTO stars are commonly used tracers owing 
to their large numbers and the ease by which they can be photomet- 
rically identified (e.g. lBell et al .120081) . The absolute magnitudes of 
such stars centre around Mr ~ 4.5 but have a wide range of values 
(o"Mr ~ 0.9 mag) which limits the accuracy to which the density 
profile can be estimated. BHB stars are superior distance estimators 
(o"i\/r ~ 0.15), but are significantly scarcer than main sequence 
stars. Moreover, they suffer from contamination by blue straggler 
(BS) stars due to their similar colours. BHB and BS stars can be 
distinguish ed by the ir Balmer line profiles (e.g. Kinman et alj 19941 ; 
lYannv et a l. 2000; Sirko et al.ll2004l; IClewleve t al. 2002), but this 
requires spectroscopic information. Whilst spectroscopic samples 
can c leanly identify BHB stars, the variety of re s ults on flattening 
(e.g. |Prestonetal.lll99ll; ISluis & Amoldl Il998l; lOe Propris et a .1 



I2OIOI) and density fall-off (e.g. lXue et al.ll2008l ; lBrown et al.ll201ol) 
suggests that the completeness biases are difficult to understand and 
control. 

An independent constraint on the density profile of the stellar 
halo is provided by t he velocity di s tributi on of the halo stars. A 
kinem atic analysis bv lCarolloetallfeOOTh (see also lCarollo et al.l 
l201(t) suggests that the stellar halo comprises of two components 
with different density profiles and metallicities. The authors find 
that the density profile becomes shallo wer beyond 15 — 20 kpc. 
This is in s t ark co ntrast to the studies bv lWatkins et al.l ( I2OO9I) and 
ISesar et al. who find a significantly steeper density profile 

beyond ~ 30 kpc from the distribution of RR Lyrae stars in SDSS 
Stripe 82. A caveat to the interpretation of kinematic studies is that 
the density distribution is not measured directly but rather modelled 
by assuming a dark matter halo potential. 

Therefore, the present state-of-play is distressingly inconclu- 
sive and a further attack on the problem of the shape of the stel- 
lar halo is warranted. In this study, we introduce a new method to 
model both BHB and BS stars based on photometric information 
alone. We make use of the SDSS DR8 photometric data release 
which has now mapped an impressive ~ 14, 000 deg^ of sky with 
both Northern and Southern coverage. In contrast to previous work, 
we combine both the wide sky coverage of the SDSS with the ac- 
curate distance estimates provided by the BHB stars to model the 
density profile of the stellar halo out to ~ 40 kpc. 

The paper is aixanged as follows. In §2.1, we describe the 
SDSS DR8 photometric data and our selection criteria for A-type 
stars. The remainder of §2 introduces the probability distribution 
for BHB and BS membership based on colour alone and outlines 
the absolute magnitude-colour relations for the two populations. In 
§4, we describe our maximum likelihood method to determine the 
density profile of the stellar halo and in §4 we present our results. 
Finally, we draw our main conclusions in §5. 



2 A-TYPE STARS IN SDSS DATA RELEASE 8 



2.1 Data Release 8 (DR8) Imaging 



The Sloan Digital Sky Survey (SDSS ; l York et al.l2000t) is an imag- 
ing and spectroscopic survey covering roughly ~ 1 /4 of the sky . 
Imaging data is obtain ed using a CCD c amera jGunn et aL I I 19981) 
on a 2.5 m telescope jGuim et aLll2006l) at Apache Point Obser- 
vatory, New Mexico. Images are obtained simultaneously in five 
broad optical bands (ugriz; Fukugita et al. 1996). The data are pro- 
cessed thro ugh pipelines to measure photometric and astrometric 
properties (Luptonetal. 2001 ; Smith et al. 20021; Stoughton et al.l 



I2OO2I ; IPier et al.ll2003l ; llvezic et alj|2004l ; iTucker et al.ll2006h . The 
SDSS DR8 release contains all of the imaging data taken by the 
SDSS imaging cam era and covers over ~ 14, 000 deg^ of sky 
jAihara et al.l201 ll) . We select objects classified as stars with clean 
photometry. The magnitudes and colours we use in the following 
sectio ns have been corrected for extinction following the prescrip- 
tion of lSchlegel etal] ( Il998h . 

In the top panel of Fig.[T| we show the sky coverage of SDSS 
data release 8 (DR8) in equatorial coordinates. For comparison, the 
sky coverage of the SDSS data release 5 (DR5) is indicated by the 
darker grey region. The more recent SDSS data release covers both 
Northern and Southern latitudes. We exclude latitudes |ti| < 30°, 
so as to concentrate on regions well away from the plane of the 
galaxy. Over the distance range probed by this work, the SDSS 
footprint (|bj > 30°) covers approximately 20% of the total vol- 
ume of the stellar halo. In this study we use blue horizontal branch 
(BHB) and blue straggler stars (BS) to map the density profile of 
the stellar halo. These A-type stars are selected by choosing stars 
in the following region in colour-colour space; 



0.9 < u-g < 1.4 
-0.25 < g-r < 0.0 



(1) 



This selection is similar to other work using A-type stars (e.g. 
lYannv et al]|2000l ;[ Sirko et al .l2004h and is chosen to exclude main 
sequence stars, white dwarfs and QSOs. The bottom left hand panel 
of Fig. [T] highlights the colour selection box. Whilst we assume 
that all of our selected stars are BHBs or BSs, there may be a non- 
negligible contribution by variable stars, such as RR Lyrae. We use 
multi-epoch stripe 82 data to estimate the fraction of variable stars 
in the same magnitude and colour range as our sample. W e use the 
light curve catalogue compiled by iBramich et al.l ( |2008|) and clas- 
sify va riable stars according to the criteria outlined in Isesar et al.l 
toot . The resulting fraction of variable stars is ~ 5%. This small 
fraction of non-BHB and non-BS stars will therefore make little 
difference to the results of this work. 

In the bottom right hand panel of Fig. [T] we show the error 
in u — g and g — r colours as a function of g band magnitude. The 
photometric errors in u — g are larger than in g — r, especially at 
fainter magnitudes. This can be compared with the typical separa- 
tion between BHBs and BS stars in u — g colour (see bottom panel 
of Figure|2]l which ranges from 0.05 to 0.1 mag. Mean photometric 
error o^u^g) reaches 0.05 at g ~ 18.5 and beyond that the value 
rapidly increases. Accordingly, in this study we only select stars in 
the magnitude range 16 < g < 18.5. This corresponds to a dis- 
tance range of ~ 4 — 40 kpc for typical absolute magnitudes of 
BHB and BS stars. Note that BHB and BS stars have different ab- 
solute magnitudes, and so span separate, but overlapping, distance 
ranges (see Section l23] l. 



© 201 1 RAS, MNRAS 000,[T]{T3] 



The Milky Way Stellar Halo 3 





Figure 1. Top panel: The SDSS DR8 footprint in equatorial coordinates. The solid and dotted lines show |ti| = 30° and 6 = 0° respectively. The darker grey 
area shows the sky coverage of the SDSS data release 5 (DR5). The more recent DR8 sample has both Northern and Southem sky coverage. Bottom left panel: 
The colour selection in u — g and g — r used to select BHB and BS stars. Our sample consists of = 20290 stai's in the magnitude range 16 < g < 18.5. 
Bottom right panel: The median (dots), 5th and 95th percentiles of the photometiic error m u — g (black) and g — r (red) colours as a function of g band 
magnitude. Beyond g ~ 18.5, the en'ors in u — g increase fairly rapidly. 



2.2 Ridgelines in Colour-Colour Space 

We seek to measure the centroids of the BHB and BS loci in colour 
space as well as their intrinsic widths. Naturally, this can only be 
done provided there exists a robust classification of the A-type 
stars according to their surface gravi t y. It has been s hown (e.g. 
IClewlev et allbool ISirko et al.ll2004 fxue et aLlllOOSl) that BHB 



and BS stars can be separated cleanly on the basis of their B aimer 
line profiles. We proceed by selecting A-type stars from the spectral 
SDSS data release 7 (DR7) catalogue within the same colour range 
as our DR8 photometric sample. Restricting the sample to high 
signal-to-noise (S/N) spectra in the magnitude range 16 < 5 < 17, 
we fit the Balmer lines H^, Hs and Hp with Sersic profile of the 
form, 



y — 1.0 — a exp - 



- 2:0 



(2) 



where xo and a give the wavelength and the line depth at the line 
centre respectively. The parameters b (the scale width) and c are 
related to the line width and line shape respectively. The relation 
between combined line widths and line shapes of the three Balmer 
lines Hj, Hg and H/s are shown in the top panel of Fig. |2] The 
BHBs (blue points) and BSs (red points) are clearly separated in 
this diagram and the decision boundary is indicated by the dashed 
black line. We use this spectral classification to pinpoint the loci of 
the two populations in the u — g, g — r colour-colour space. In our 
analysis, these "ridgelines", i.e. approximate centres in u — p as a 



function of g — r, are defined by third order polynomials: 

(w-ff)BHB = 1.167 - 0.775((?-r)- 1.934(ff-r)2 
+9.936{g-rf, 
(w-5)bs = 1.078 - 0.489((7-r)+0.556(ff-r)" 
+ 13AU{g-rf, 



(3) 



for —0.25 < .g—r < 0.0. The ridgelines are shown by the thick blue 
and red lines in the bottom panel of Fig.|2] The green line indicates 
the approximate boundaiy between BHBs and BSs in u — g, g — 
r space. In addition, we calculate the intrinsic spread of the two 
populations about their ridgelines. We find o"bhb,o(m — <?) = 0.04 
and cFBSfiiu- — g) = 0.045. Fig. |2] makes it apparent that, even 
for brightest stars, photometric information alone is not enough to 
separate BHB and BS stars. Therefore, for a given star we define the 
probability of BHB or BS class membership based on its distance 
in colour-colour space from the appropriate locus. We assume that 
both populations are distributed in a Gaussian manner about their 
ridgelines. We model the conditional probability of measuring u—g 
and g — r colours, given the star of each species, as 



p(ugr I BHB) oc exp 



[{u-g) - {u-gfl^^Y 



^"bhb 



p{ugr I BS) oc exp I — 



[{u -g)-{n-g)%sy 



(4) 
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Figure 2. Top panel: The scale width, 6(7(5/3), and line shape, c(7(5/3), of 
the Balmer lines Hj, Hs and i?^ for bright A-type stars taken from the 
SDSS data release 7 (DR7) spectral catalogue. These stars are selected in 
the same colour range as our DR8 photometric sample and have magnitudes 
in the range 16 < g < 17. The blue and red points denote BHB and BS 
stars respectively. The black dashed line shows the apparent separation of 
these stars based on the Sersic line profiles of the Balmer series. Bottom 
panel: The 'ridgelines' for the BHB and BS stars in the colour space (u — 
g, g—r). The thick blue and red fines sfiow tfie third order polynomial fits to 
these loci in colour space. The green line indicates the approximate border 
between the two populations. 



Note that, in fact, these probabilities are also functions of g — r 
since the centre of the Gaussian distribution, (it — g)" varies with 
the g — r colour (see eqn. [3). The dispersion about the ridgehne 
centre depends on the intrinsic width and the photometric errors in 
u~g, a = ,/ + • The colour-based posterior probabilities 



("-9)' 

of class membership are then 



P(BHB I ugr) 
P(BS I ugr) 



p{ugr I BHB) iVBHB 



p(ugr I BHB) iVsHB + p{ugr \ BS) A^bs 

pjugr I BS) Ai'BS 

p(ugr I BHB) iVBHB + p{ugr \ BS) A''bs 



(5) 



The total numbers of stars A'^bhb and TVbs in a given colour range 
can then be found iteratively by integrating equations (|5)- In Ta- 
ble [T] we give the fraction of BHB and BS stars in five g — r 
colour bins. The fraction ranges from /bhb ~ 0.15 at redder 
colours to /bhb ~ 0.6 at bluer colours. This is in good agreement 
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Figure 3. The distribution of colours in (u — g)* space. Each row shows 
the distribution for a different range in g — r colour A two component 
model is fit to the data. The overall model is shown by the black line. Indi- 
vidual Gaussians are shown by the red (BS) and blue (BHB) lines respec- 
tively. Stars with [u-g)* < are dominated by BS stars whilst stars with 
{u — g)* > are dominated by BHB stars. The ratio between the amplitude 
of the Gaussians gives an estimate of the overall number ratio between the 
two populations. 



wit h the overa l l BHB to BS ratios estimated bv lBell et all j2010t) 
and lXue et all ( |2008|) who use similar magnitude ranges. Figure |3] 
demonstrates the evolution of the separation between the two pop- 
ulations in colour space. For illustration purposes, the SDSS u — g 
colour is transformed into (u — g)* using the following relation: 



{u-g)* = (u-g) - {t- 



/border 



(6) 



Here, {u-g)tri,, = 1.223 



0.632(,g-r) - 0.689(5 -r)^ + 



11.690(,g-^ ) is defined by the approximate boundaiy line between 
BHB and BS stars shown by the green line in Fig. |2] In these new 
coordinates, the curved shape of the decision boundary becomes a 
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/bhb 


/bs 


-0.05 < g~ 


-r < 0.00 


5189 


0.154 


0.846 


-0.10 < g~ 


r < -0.05 


4973 


0.231 


0.769 


-0.15 < g~ 


r < -0.10 


4151 


0.346 


0.654 


-0.20 < g- 


r < -0.15 


3564 


0.536 


0.464 


-0.25 < g- 


r < -0.20 


2413 


0.613 


0.387 



Table 1. The fraction of BHB and BS stars in different colour bins. We give 
the g — r colour range, the total number of stars, the estimated fraction of 
BHB stars and the estimated fraction of BS stars. 



straight line. In Fig.|3] we fit two Gaussian distributions to the dis- 
tribution in (u — g)* for five bins in g-^. This two component model 
fits the overall distribution very well. The ratio between the ampli- 
tudes of the two Gaussians varies with g—r colour. The fraction of 
BHB stars increases towards bluer colours, in good agreement with 
our estimates in Table [T] 



2.3 Absolute Magnitudes 

Let us now derive a relationship between the absolute magni- 
tude and colour of BHB and BS stars. BHB stars are intrinsically 
brighter and their absolute magnitude varies little as a function of 
temperature (colour) or metallicity. In comparison, BS are intrinsi- 
cally fainter and span a much wider range in absolute magnitude. 

The absolute magnitudes of BHB stars are calib rated us- 
ing st ar clusters with SDSS photometry published by I An et al.l 
( 120080 . Ten star clusters have prominent BHB sequences; M2, 
M3, M5, M13, M53, M92, NGC2419, NGC4147, NGC5053 and 
NGC5464E 

The density distribution of absolute magnitudes of 
BHB stars in these clusters is shown as a function of g—r colour in 
the left hand panel of Fig.|4] There is little variation of the absolute 
magnitude with colour for BHB stars. Afg changes by 0.2 (from 
0.65 to 0.45) in the —0.25 < g — r < range. The inset panel 
shows the variation of absolute magnitude about the Mg(BHB) vs. 
g—r trend (Mg ). The spread of this distribution is ~ 0.1, indicating 
that there is a tight relation describing the BHB absolute magnitude. 
The star clusters have metallicities typical of halo stars and ranging 
from —2.3 < [Fe/H] < —1.3, but we find no obvious trend with 
metallicity. 

BS stars are less common in globular clusters than BHB stars. 
Instead, to calibrate the absolute magnitudes of BS stars we make 
use of stars in Stripe 82 belonging to the Sagittarius stream. The 
distance to the stre am in the right ascensi on range 25° < a° < 40° 
was estimated by IWatkins et al.l ( |2009|) using RR Lyrae stars as 
Dsgi ~ 26.1 ± 5.6 kpc. In the middle panel of Fig.|4l we show the 
absolute magnitude of stars in Stripe 82 between 25° < a° < 40° 
as a function of colour. The density contours are constructed by us- 
ing stars outside of the range in right ascension as a background 
and computing the density contrast. An obvious plume of BS stars 
is apparent in the density plot. This is shown by the contour levels 
extending off the main sequence. The solid red line shows the es- 
timated absolute magnitude colour relation for these BS stars. For 



1 We adopt distance moduli for NGC2419 and NGC4147 of 19.8 and 16.32 
respec tively. These differ form the values given in Table 1 of lAnetalJ 
120081) . We find that these revised values are in better agreement with the 
colour magnitude diagrams of the clusters. 



comparison, we show the abs olute magnitude vers us colour rela- 
tion for BS stars estimated bv lKinmanetalJ l ll994l) . This relation 
is converted from Johnson-Cousins photometry (UB V) to Sloan 
photom etry (ugr) using the transformation derived in Ijester et al.l 
( l2005i) . The different coloured dots show the relation for differ- 
ent metallicity stars. The absolute magnitude calibration derived 
for the BS stars in the Sagittarius stream is almost identical to 
the iKinman et al ] < 19941) relation for BS stars with metallicity 
[Fe/H] = —1.5. This is in good agreement w ith the metallicity of 
the stream stars found bv lWatkins et al.l ( l2009l) ([Fe/H] = -1.43). 

We compute the spread of absolute magnitudes for each colour 
bin as ~ 0.5. This dispersion takes into account the distance 
errors (Dsgr = 26.1 ± 5.6 kpc) and encompasses a range of metal- 
licities (see dashed red lines in Fig. |4j. Hence, we conclude that 
our calibration for BS absolute magnitudes does not have a strong 
metallicity bias. Note that the 'average' BS absolute magnitude is 
~ 2.5, approximately 2 magnitudes fainter than the BHB stars. 
This is entirel y consistent with th e abs olute magnitudes of h alo BS 
stars found bv lYannvet^ ( |2000|) and lClewlev et al.l ( |2004|) . 

The resulting absolute magnitudes of BHB and BS stars as a 
function of g — r colour are; 

Mg(BHB) = 0.434 - 0.169(,g-r-) + 2.319(5-r)' 
+20.449(5-r)'^ + 94.517(5-r)'', 
Mg(BS) = 3.108 + 5.495(3-r), (7) 

which are valid over the colour range —0.25 < g — r < 0.0. For 
BHB stars, this allows us to estimate accurately the distances of 
BHB candidates. For BS stars, the relation is only approximate and 
the scatter for each colour interval needs to be taken into account 
for all distance estimates. 

In the right hand panel of Fig.|4] we show the estimated radial 
distributions for the BHB and BS populations. We select high prob- 
ability BHB and BS stars by using the membership probabilities 
defined in eqn. (|5]l. 'High' probability BHB/BS stars are defined as 
those for which we are 68% (or la) confident of BHB/BS mem- 
bership. The BHB stars probe a radial range from r ~ 10 kpc to 
r ~ 45 kpc. BS stars, which have fainter absolute magnitudes, only 
probe out to ~ 30 kpc. However, there is a large degree of overlap 
between the two populations between 10 and 30 kpc. 



3 MAXIMUM LIKELIHOOD METHOD 

In this section, we outline the maximum likelihood method used to 
constrain the density profile of the stellar halo. An important as- 
sumptions in the modelling is that the BHB and BS stars follow the 
same density distribution, modulo an overall scaling. The number 
of BHB stars and BS stars in a given increment of magnitude and 
area on the sky is then described by 

ANsHBimg - Mf™ , £, b) = pbhbP('t1s - A/f ™, £, 6)i3iHB 

X ilnlO Am„cosfeA.^Ab 
5 

AiVBs(mg-A//f ,A6) =PBsp(mg-A//f ,^,6)i3|s (8) 

X ilnlO Am„cos6A^A6. 
5 

Here, the distance increment AD has been converted into the ap- 
parent magnitude increment via the relation AD = ilnlO DArrig. 
The normalising factors, Pbhb = A'bhb/Vbhb and Pbs = 
Nbs /Vbs are found by performing volume integrals over the SDSS 
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Figure 4. Left Panel: The colour-absolute magnitude relation for BHB stars derived from star clusters published in lAn et allfcOOSh . A polynomial of order 4 is 
fit to the BHB stars in eleven star clusters: M2, M3, M5, M13, M53, M92, NGC2419, NGC4147, NGC5053 and NGC5466 (solid blue line). The grey contours 
indicate the density of stars within the colour-absolute magnitude region (thicker lines indicate higher densities). The ridgelines for individual clusters are 
shown by the red/black lines. These are colour coded by metallicity: red to black goes from more metal rich ([Fe/H] -1.3) to more metal poor ([Fe/H] -2.3). 
The inset panel shows the distribution of values around the derived relation indicating a small degree of scatter Middle panel: The density of BS stars in Stripe 
82 belonging to the Sagittarius stream. The solid and dashed red lines give the absolute magnitude colour relation and the estimated dispersion ipM^ ~ 0.5). 
The distan ce to this portion of t he Sagittarius stream is estimated from Watkins et al. (2003) as Dggr = 26.1 it 5.6 kpc. The coloured dots show the relation 
derived bv lKinmanetalJ<1994b for different metallicity BS stars in star clusters. Right panel: The radial distribution of the BHB and BS star populations. The 
blue fine shows the distribution for high probability BHB stars (P(BHB) > 0.7). The red shaded region shows the distribution for high probability BS stars 
(P(BS) > 0.7) where the uncertainty in the absolute magnitudes is taken into account. 



DR8 sky coverage and over the required magnitude range 



vbs(m; 



j j j [pim,-Mf™J,b)DlnB 



X -lnlOdm„ cosbd^dfe 
5 



BS 



X -InlOdmo cos&d^d& 



(9) 



These normalising integrals depend on the absolute magnitude of 
the stars and hence on the g — r colour (Mg = Mg{g — r) from 
eqn. |7j. Note that the values of Vbhb and Vbs play a minor role 
in identifying the maximum likelihood model given the choice of 
parameters, but are important when evaluating the performances of 
different model families. We also assume that our sample consists 
only of BHB and BS stars so the total number of stars is the sum 
of these two populations, A''tot = A'^bhb + A^'bs where A^bhb ~ 
/BHB-/Vtot and A'bs = /es-^tot- The overall fraction of BHB to 
BS stars varies as a function of g—r colour, as shown in the previous 
section. 

Combining equations l[4j and l[8j gives the number of stars in 
a cell of colour, magnitude, longitude and latitude space 

AN = p{ugr I BHB)A7Vbhb +p(M.gr | BS)AiVBS 

= Ntotviugr, I, fe) A(3-r) Am„ AMfocos&ilnlO, (10) 

5 

where the stellar density is 

v{ugr,l,h)=p{ugr \ BHB) ^plm^ -Aff™, ^, 6)73|hb + 



p{ugr I BS)^p(m,-Aff ,^,6)i3|s 



(11) 



Here, each star is assigned a 'BHB distance' (Dbhb) as well as a 
'BS distance' (Dbs)- The colour probability functions weight the 
contribution of each star to the BHB density or the BS density. For 
simplicity, we group all the stars into five bins in g—r of width 0.05 





a 


1 


ln(£) X 10'' 


(j/tot 


with V&S? 


20290 
15403 


^■O'-'-0.05 
2.90t^;f 


n cc-l-0.02 


-17.1171 
-12.2917 


0.38 ±0.01 
0.22 ± 0.02 


yes 
no 



Table 2. A summary of our best-fit oblate power-law models with and with- 
out the Virgo Overdensity and the Sagittaiius stream. We give the total num- 
ber of stars, the model parameters, the average log-likelihood value for the 
model and cr/tot. 



mags. Thus, stars in each g — r bin have the same normalisations 
(Vbhb, Vbs), fraction of BHB/BS stars (/bhb, ,/bs) and absolute 
magnitudes (M™^ , Mf^). We take into account the uncertainty 
in the BS absolute magnitude by convolving the BS number density 
with a Gaussian magnitude distribution. This distribution is centred 
on the estimated absolute magnitude {Mf^ = Mf^ (g—r)) and has 
a dispersion of — 0.5 (see Fig.|4j. 

The log-likelihood function can then be constructed from the 
density distribution. 



log£ = log \u{rng,ugr \ t ,b'') cosb\ 
1=1 



(12) 



The number of free parameters constrained by the likelihood func- 
tion depends on the complexity of the model stellar-halo density 
profile. The log-likelihood is maximised to find the best-fit param- 
eters using a brute-force grid search. 



4 RESULTS 

In this section, we outline the results of applying our maximum 
likelihood method to the sample of A-type stars selected from the 
SDSS DR8. We consider in turn a number of simple density pro- 
files with constant flattening - single power-law, broken power-law 
and Einasto - before examining the case for refinements, such as 
triaxiality, radial variations in shape and substructure. 
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Figure 5. Left panels: Density plots in equatorial coordinates. Right panels: Density plots in Galactic coordinates. The top and middle panels show the density 
plots for the the data and single power-law model respectively. The bottom panels show the residuals of the best-fit single power-law model. The dark red 
regions show the obvious overdensities of Virgo and Sagittarius. The dashed lines indicate the regions of sky removed in the maximum likelihood procedure. 
The Virgo Overdensity and parts of the leading tail of the Sagittarius stream are apparent in the North Galactic Cap whereas another portion of the Sagittaiius 
stream is found in the South Galactic Cap. Note that these have been excised when calculating the best-fit single power-law model. The regions away from 
these known overdensities are reasonably well fit by a smooth, power-law density model. 
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Figure 6. The maximum hkelihood contours for the flattening q and power- 
law index a of our single power-law halo models. The black (red) lines 
show the 1(T, 2(7 and 3(t contours when stars in the region of the Virgo 
overdensity and Sagittarius stream have been included (excluded). The inset 
panel illustrates that the difference in maximum likelihood parameters is 
relatively small whether the overdensities are included or excluded. 



4.1 Single Power-Law Profile 

First, let us consider a simple power-law density model of the form, 

p(r,)cxr-", rl = + + z\-^ (13) 

The parameters q and a describe the halo flattening and the power- 
law fall-off in stellar density respectively. Oblate density distribu- 
tions have q < 1, spherical q = 1 and prolate g > 1. 

We summarise the best-fitting single power-law model param- 
eters in Table |2] Using the entire sample, we find the maximum 
likelihood model parameters of a = 2.6 and q = 0.65. We repeat 
the analysis by excising stars in the regions of the Virgo Overden- 
sity and Sagittarius stream, which amounts to removing stars in the 
regions defined by 

< X < 30, X = 63.63961^2(1 -sinfe). (14) 

This is the mask introduced by ISell et all ( |2008|) to remove stars 
belonging to the Virgo overdensity (as well as parts of the leading 
tail of the Sagittarius stream). Another portion of the Sagittarius 
stream is located in the Southern part of the sky (see Fig. |5) and is 
removed by 

0°<a<50°, -30°<(5<0°. (15) 

By discarding these large overdensities, we find a slightly steeper 
power law and a more flattened shape with a = 2.9 and q — 0.53. 
There is also a substantial increase in the maximum likelihood, im- 
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Figure 7. Likelihood contours for the broken power-law models. Left panel: The contours for the inner (Oin) and outer (Oout) power laws. The fall off steepens 
at larger radii. Right panel: The contours for the break radius (rj,) and flattening (q). The maximum likeUhood solution favours a flattening of g = 0.59 and a 
break radius of rj, = 27 kpc with imier and outer power-laws of din = 2.3 and Qout = 4.6 respectively. 



plying that the fit is indeed affected by the presence of these large 
overdensities. In Fig |6l we show the maximum likelihood contours 
for the model parameters q and a for both cases. The contours en- 
compass the la, 2a and 3a regions respectively. We find that our 
likelihood function has a well defined peak. The maximum likeli- 
hood model parameters (a = 2.9, q = 0.53 excluding overdensi- 
ties) are in good agreement with some of the previous work assum- 
ing a single power-law model for the stellar halo (e.g. lYannv et al.l 
l200d : iNewberg & Yannvll2006l : | Juric et al ]|2008h . This exercise in- 
dicates that the two large overdensities induce a bias in deduced 
model parameters and, therefore, in the subsequent sections, we re- 
move stars in the vicinity of the Virgo Overdensity and Sagittarius 
stream in all our calculations. 

We use our maximum likelihood smooth oblate halo model 
(a = 2.9, q = 0.53) to show the residuals of the data minus 
model on the sky. In Fig. [5] we compare the data and model in 
both equatorial and galactic coordinates. The top and middle pan- 
els show the data and model on the sky (the magnitude distribution 
has been collapsed) whilst the bottom panels show the residuals of 
the model. The Virgo overdensity is the most obvious feature lo- 
cated at (q ~ 190°, S ~ 30°). This overdensity covers a substan- 
tial fraction of the North Galactic cap. A portion of the Sagittarius 
stream can be seen at (a ~ 20°, 5 ~ —25°) in the South Galac- 
tic cap. The residuals for these features reach up to approximately 
three times the model values. We can estimate the total fraction 
of stars residing in these overdensities from the excess numbers of 
stars in these regions of the sky. Approximately ~ 5.5% of our total 
sample (A'tot = 20290) reside in the Northern Virgo-l-Sagittarius 
overdensity whilst less than 1% occupy the Southern portion of the 
Sagittarius stream. The fraction of stars in these overdensities is rel- 
atively small as they are located in the vicinity of the Northern and 
Southern Galactic caps where the density of stars is small. How- 
ever, as the relative difference (i.e. (data-model)/model) between 
the data and model is large in these regions, these overdensities can 
influence the likelihood values. 



4.2 Broken Power-Law Profile 

We relax our models to allow a change in the steepness of the den- 
sity fall off and consider broken power-law models of the form. 



In Fig. [7] we show the maximum likelihood contours for the 
inner and outer power-laws (left hand panel) and break radii and 
flattening (right hand panel). A model with break radius rb — 27 
kpc is preferred with slopes of ain — 2.3, Qout ~ 4.6 respectively. 
A steeper power-law is favoured at larger radii while the power-law 
within the break radius is shallower. Note that the break radius is an 
ellipsoidal distance and only corresponds to a radial Galactocentric 
distance in the Galactic plane (z — 0). The model is slightly less 
flattened (g ~ 0.6) than the single power-law model but, even with 
an additional two parameters, an oblate halo model is favoured. 

The maximum likelihood value for a broken power-law 
model is significantly larger than for a single power-law model 
(-21n(£sPL/£BPL) ~ 400) (see Tabled. Our broken power- 
law model is in very good agreement agreement with the results of 
IWatkins et al. U2009I) . who use a sample of RR Lyrae stars to probe 
the density distributio n of the stellar halo out to '-^ 100 kpc (see 
also lSesare"tal]|20I0l) . 



4.3 Einasto Profile 



p(r,) oc 



rq > Tb. 



(16) 



The Einasto profile jEinasto & Haujl989h is often u sed to describe 
the density distribution of dark matter haloe s (e.g. iGraham et al.l 
I2OO6; Merritt et al. 2006: Navarro et al.lboTol) . The Einasto model 
is given by the equation 

In [pirq)/p{r,a)] = -d„[(r,/reff)^/" - 1]. (17) 

The shape of the density profile is described by the parameter 
n. Density distributions with steeper inner profiles and shallower 
outer profiles are generated by large values of n. For example, dark 
matter haloes with 'cuspy' inner profiles typically have values of 
n ~ 6. The parameter rf„ is a function of n. For n ^ 0.5 a good ap- 
proxi mation is given by d„ = 37i— 1/3 + 0. 0079/n JGraham et al.l 
I2OO6I) . 

This profile allows for a non-constant fall-off without the need 
for imposing a discontinuous break radii. We find the maximum 
likelihood solutions for q, n and r-cfi. Our best-fit Einasto model 
has parameters q = 0.58, n = 1.7 and r^fi — 20 kpc. The slope 
of the density profile varies rapidly with radius as indicated by the 
relatively small value of n. 
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Figure 8. The magnitude distribution of our DR8 data sample. The black points give the distribution of the data where the en'or bars are Poissonian. Top-left 
panel: Single power-law models. The solid and dashed red lines show the best fit oblate and spherical models respectively. Error bars incorporate Poissonian 
uncertainties and the spread of absolute magnitudes for BS stars. A flattened model provides a much better fit to the data. Top-right panel: The solid red, blue 
and green lines show the distributions for the best-fit single power-law, broken power-law and Einasto models respectively. The latter two provide a better 
representation of the data. Bottom-left panel: The magnitude distribution for the most probable BHB stars (with Pbhb > 0.7). Bottom-right panel: The 
magnitude distribution for the most probable BS stars (with Pbs > 0.7). 



4.4 Model Comparisons 

We now test how accurately our models represent the observed 
magnitude distribution of the data. Using Monte Carlo methods, 
we create a distance distribution according to the model density 
profile. The fake data is given a g — r colour distribution drawn 
from the real data, which can then be converted into absolute mag- 
nitudes (see eqn.|7]l. The ratio of BHB and BS stars in each colour 
bin is chosen to match the values given in Table [T] In the case of 
the BS stars, absolute magnitudes are determined from the g~r 
colour by drawing randomly from a Gaussian distribution centered 
on the estimated value (from eqn.|7} with a dispersion of 0.5 mags. 
The resulting (apparent) magnitude distributions for our models are 
compared with the data in Fig. [8] 

In the top-left hand panel, we show the magnitude distribution 
for our single power-law oblate model with the solid red line. The 
error bars take into account the Poisson uncertainty as well as the 
uncertainty spread of absolute magnitudes for BS stars. The dis- 
tribution of magnitudes for our DR8 A-type star sample is shown 
by the black points. There is reasonably good agreement but there 
is a notable deviation at fainter magnitudes. For comparison, we 
show the maximum likelihood spherical model by the red dashed 
line (with a ~ 2.7), which provides a very poor fit to the data. 
The top-right hand panel shows the magnitude distributions for 
the single power-law, broken power-law and Einasto profiles by 
the solid red, blue and green lines respectively. The broken power- 
law and Einasto models provide a better representation of the data 



than a single power-law, although there are still discrepancies at the 
faintest magnitudes. 

In the bottom panels of Fig. [8] we show the magnitude dis- 
tributions for the most probable BHB and BS stars. We select stars 
with Pbhb > 0.7 as 'BHB' stars and stars with Pbs > 0.7 as 'BS' 
stars (see section 12. 3t . Of course, this is for illustration purposes 
only, as a clean separation of BHB and BS stars requires spectro- 
scopic classification. While the BS stars are adequately described 
by a single power-law model, this is a poor description of the BHB 
stars. This is unsurprising, as the BS stars cover a smaller distance 
range and barely populate distances beyond the break radius. The 
bottom-left panel clearly shows the need for a more steeply declin- 
ing density law at larger radii. 

In Fig. [21 we show the distribution of (probable) BHB stars 
in spherical shells with the solid black line. Here, we only consider 
the most likely BHB stars (with Pbhb > 0.7) as for these we can 
accurately estimate their distance. The red, blue and green lines 
show the radial distributions for our best fit single power-law, bro- 
ken power-law and Einasto models respectively. The single power- 
law is a poorer description of the data, whilst the broken power-law 
and Einasto models both provide better representations of the radial 
distribution. 



The Bayesian evidence is the integral of the likelihood values over the 
parameter space (assuming a uniform prior). E j C{6)A9, where d is 
the model parameter vector. 
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Model Pai-ameters Np ln(£) X 10* -21n{£max//;) ln(_B/£;max) o-/tot (w/o V & S) fT/tot(V&S) 



SPL- 


spherical 


a 


2 7+0.05 _ 
^- ' -0.05- 1 — 


1.0 


1 


-12.5516 


5606 


-2801 


0.44 ±0.01 


0.49 ±0.01 


SPL- 


oblate 


■ 


2 9+0.04 _ 
^•^-0.06' 1 ~ 


'^•^■'-o.oi 


2 


-12.2917 


408 


-206 


0.22 ±0.02 


0.38 ±0.01 


SPL- 


triaxial 


■ 


2 n+0.05 _ 
^•='-0.05' 1 ~ 


n 50+0 02 „ 7, +0.03 


3 


-12.2818 


210 


-110 


0.20 ±0.02 


0.34 ±0.01 


SPL- 


q = q{r) 


■ = 


2 0+0.05 

^•='-0.05' y ~ 


0.531°;°^, ro > lO^kpc 


3 


-12.2917 


368 


-206 


0.22 ±0.02 


0.38 ±0.01 


BPL- 


oblate 


rb = 


27tikpc, g = 




4 


-12.2713 








0.21 ±0.02 


0.36 ±0.01 



ctin = 2.3+^;J,aout = 4.6I^;| 

Einasto - oblate n = 1.7^0,1, r^ft = 20t{'!^^pc, q = 0.58^0 gl ^ -12.2757 88 -45 0.22 ± 0.01 0.37 ±0.01 



Table 3. A summary of our best-fit models. We give the type of model, the best-fit parameters of the model, the number of free parameters, the maximum 
log-likelihood, the difference in likelihood relative to the maximum likelihood model, the log evidence^ ratio relative to the maximum likelihood model and 
cr/tot both with and without the Virgo and Sagittaiius overdensities. Parameters which are kept fixed are highlighted in bold. 



4.5 Refinements: Triaxiality and a Radially Dependent Shape 

A natural question to ask is whether further refinements might pro- 
vide a still more accurate description of the data. Up to now, we 
have assumed spheroidal halo models with a constant flattening 
with radius. 

First, we consider whether triaxiality makes any improvement. 



(triaxial-oblate)/oblate 



2 I 2 - 
— X + y p 



2 , 2 - 
+ z q 



(18) 



We fit single power-law triaxial models to the data and obtain max- 
imum likelihood parameters of q = 0.50 and p — 0.71. The likeli- 
hood value increases relative to an oblate model with the inclusion 
of this extra parameter ( — 21n(£obiato/Ariaxiai) ~ 200). The mag- 
nitude distribution of the triaxial model is largely the same as the 
oblate model. We inspect the difference between these two mod- 
els in equatorial coordinates in Fig. [TO] The dot-dashed lines in- 
dicates the regions of sky with low galactic latitudes \b\ < 40°. 
Regions with latitudes below |6| < 30° have been removed. The 
triaxial model is notably overdense relative to the oblate model 
in the regions centered on (a, 5) = (125°, 50°) and (a, 5) = 
(60°, 10°). The latter region, located in the Southern part of the 
sky, is close to the portion of the Sagittarius stream excised in our 
best fit model. The f ormer is coincident with the Monocerus ring 
jNewberg et alj2002h. wh ich ca n be identified, for example, in Fig- 
ure 1 of Belokurov et al. liooi). We suggest that the apparent tri- 
axiality may well be caused by the presence of these overdensities. 
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Figure 10. The residuals of our maximum likelihood triaxial and oblate 
models. The dot-dashed region indicates the boundary between low < 
40°) and high (|fe| < 40°) Galactic latitudes. Note the difference in scale 
that is used for this figure to that used for Fig.|5] 



The increased flexibility of the model can cause it to 'fit' to such 
substructure and thus the increase in likelihood may be an artifact. 

Some earlier investigations have found evidence that the 
shape of the stellar halo changes from a flattened distribution at 
small er radii to an almo st spherical distribution at larger radii 
(e.g. IPreston et al. I I1991I) . We can test this claim by allowing 
the flattening n to vary with radius. Following the reasoning of 
ISluis & Arnold ( Il998i) . we make the following substitution 



+ 4 



q2j.: 



(19) 
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Figure 11. A side view of our maximum likeliliood models. Greyscale sliows density of stars in plane of Galactocentiic (x, z) at y = in four maximum- 
likelihood models. The blue, green, yellow and red contours show density levels corresponding to the 50th, 90th, 95th and 99th percentiles of the spherical 
model. The white dot marks the location of the Sun. 
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Figure 12. The scale dependence of substructure. We show the cr/tot fraction as a function of superpixel size. Individual panels show the relation for our 
various maximum likelihood models. The different line styles indicate different magnitude bins. The black lines show the relation for all the stars. The red 
lines show the relation when stars in the vicinity of the Virgo overdensity and Sagittarius have been removed. 



The halo flattening is still q at small radii but tends to sphericity 
at large radii. The scale radius ro determines the radial range over 
which this change occurs. For example, for large values of ro (e.g. 
larger than the most distant stars) the flattening is approximately 
constant over the applicable radial range. We fit single power-law 
models with a varying flattening to the data and find that large 
scale radii are preferred (ro > 10^ kpc) with an inner flattening 
of g = 0.53. This indicates that the flattening is approximately 
constant over the radial range of the d ata (out to ~ 40 kpc). This 
is in agreement w ith the deductions o f ISluis & Amoidl ( ll998l) and 
ISesar et al] ( 1201 1[) . who also found no real evidence for a varying 
shape with radius. 

In summaiy, we find that the data are well described by an 
oblate density distribution with a constant flattening of g ~ 0.6 
and a more steeply declining profile at larger radii. We summarise 
our maximum likelihood models by giving a 'side view' of the pro- 
files in Fig. [TT] A spherical model is strongly disfavoured whereas 
oblate broken power-law and Einasto models provide a good repre- 
sentation. 

Our best-fit density distribution can be used to estimate the 
total steUar mass. We find the total number of BHB stars by inte- 
grating the BHB density profile over all space in the distance range 

1 — 40kpc. The number of BHB stars can b e converted in t o a to- 
tal Luminosity using the relation derived in iDeason et al. I( l201lh 
using globular clusters: Nbus/L ~ 10"'^. Assuming a mass-to- 
light ratio of M/ L ~ 1 — 5, the stellar mass is approximately 

2 - 10 X lO^Mr:^. ISeil et al.l j200^ calculate a total stellar mass 
of ~ 3.7 X 10* Mq using main sequence stars, in good agreement 
with our estimate. Not e that the total stellar halo mass is believed 
to be ~ lO'^M© (e.g. iMorrison et all 1 1993b so we are probing a 
significant fraction of the stellar halo (~ 20 — 100%). 



4.6 The Amount of Substructure 

A rough idea of the amount of substructure present in the data can 
be attained by computing the rms deviation of the mo dels about 
the da ta (cr/tot) as defined by equations (2) and (3) in lBell et al.l 
( l2008l) : 

{cj^) = -y^{Dk-Mkf --Y^ (mI -mX , (20) 



tot (l/n)E.^fc' 
Here, Dk is the number of stars in bin fc, Mh is the expected num- 
ber from the model, is a Poisson random deviate from the 
model value and n is the total number of bins (in /, 6 and mg). 
However, with a total of ~ 20, 000 stars we can not afford to use 
a fixed bin size over the entire sky: pixels must be simultaneously 
large enough to contain ample stars and small enough to adequately 
sample the substructure. T o circumvent this problem, we use the 
Voronoi binning method o f ICappellari & CopinI (l2003h to partition 
the data in pixels on the sky. This adaptive binning method groups 
pixels together, forming 'superpixels', with the objective of obtain- 
ing a constant signal-to-noise ratio per bin. We choose the signal- 
to-noise ratio of SjN ~ 4 (assuming Poisson noise) to ensure 
the mean number of stars in each 2D bin is > 10. We also check 
that our results are not affected by the choice of signal-to-noise. 
The stars are split into three g band magnitude ranges, binned into 
~ 65, 000 1° X 1° pixels for full sky coverage which are then com- 
bined into ~ 750 superpixels using the procedure described above 
to estimate the overall cr/tot values. These are given in the last two 
columns of Table[3] As expected, our a /tot fraction is significantly 
reduced when the stars in the region of the Virgo Overdensity and 
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Sagittarius stream are removed: tlie estimated fraction of substruc- 
ture reduces from 40% to 20%. However, furtiier refinement of our 
model makes little difference to the a / tot fraction, even if the lilce- 
lihood is increased. This is indicative of the resilience of our maxi- 
mum likelihood method to relatively minor amounts of substructure 
present in the data. 

We illustrate the dependence of cr/tot on spatial scale in Fig. 
I12lbv grouping together superpixels with similar spatial scales and 
computing a/tot for each group. This shows how the fraction of 
substructure depends on the superpixel size (in deg^) for differ- 
ent maximum likelihood models, cr/tot is computed both with 
(black lines) and without (red lines) the stars in the vicinity of the 
Virgo Overdensity and the Sagittarius stream. The spherical model 
has relatively large values of a/tot over all scales, illustrating the 
poor fit of this model to the data. The single and broken power- 
law oblate models show a scale dependence, whereby larger scales 
(larger number of pixels) have an increased a/tot fraction rela- 
tive to smaller scales. Similar behaviour is also true for the triax- 
ial model but the a /tot fractions are slightly lower. As alluded to 
earlier, we suspect that this may be caused by the triaxial model 
describing some of the large-scale substructure (namely the Mono- 
cerus ring or parts of the Sagittarius stream which have not been 
excised). 

Different lines in Fig. [T2] show the dependence of the a /tot 
measure on the apparent magnitude of the tracers and hence corre- 
spond to different distances probed. The slight rise of the rms with 
magnitude indicates that higher fractions of substructure may be 
present at larger distances. This can be simply explained with the 
increase of the substructure lifetime with radius as governed by the 
orbital period of the infalling Galactic fragments. 

When interpreting these plots, it is important to bear in mind 
that the smallest superpixels are biased towards the densest parts of 
the halo, while the emptiest regions are sampled by the largest ag- 
glomerations of 2D bins. This bias could potentially explain at least 
some of the trend with size. On the scales smaller than several tens 
of degrees, cr/tot in our models is in the range 0.05 < cr/tot < 
0.2 indicating that the halo is not dominated by substructure and is 
relatively smooth. While there is evidence that the fraction of sub- 
structure increases with scale, this apparent increase in cr/tot could 
also be caused by spurious pixels. In sparse regions, the grouping 
of many pixels, each containing very few stars, poses a limitation 
to this method. 



5 CONCLUSIONS 

We have developed a new method to simultaneously model the den- 
sity profile of both blue horizontal branch (BHB) and blue straggler 
(BS) stars based on their broad-band photometry alone. The prob- 
ability of BHB or BS membership is defined by the locus of these 
stars in u — g, g — r space. We use these colour-based weights to 
construct the overall probability function for the density of the two 
populations. When applied to the A-type stars selected from the 
Sloan Digital Sky Survey (SDSS) Data Release 8 (DR8) photomet- 
ric catalog, the best-fit stellar halo models are identified by applying 
a maximum likelihood algorithm to the data. 

Based on the data from regions with 30° < |&| < 70°, we 
find that the stellar halo is not spherical in shape, but flattened. 
Spherical models cannot reproduce the distribution of the A-type 
stars in our sample and are discrepant with the data at both bright 
and faint magnitudes. Our best-fit models suggest that the stellar 
halo is oblate with a flattening (or minor axis to major axis ratio) 



of g = 0.59 with a typical uncertainty of cr, ~ 0.1. As a simple 
representation of the stellar halo, we advocate the use of a broken 
power-law model with an inner slope ain = 2.3 and an outer slope 
Qout = 4.6, the break radius occurring at about 27 kpc. This gives 
the formula 



p(r,) oc 



27 kpc, 
> 27 kpc, 



(22) 



with r-q = R + z /0.59 . For those who prefer Einasto models, 
an equally good law for the stellar halo density is 



p{rg) oc exp [-4.77[(r-5/20kpc) 



1/1.7 



-11 



(23) 



These two formulae are fundamental results of our paper, and can 
be summarised as the Milky Way stellar halo is squashed and bro- 
ken. The results are qualitatively similar to those obtained by Sesar 
et al. (2011) using main-sequence stars. There is no evidence in 
the data for variation of the flattening with radius. There is some 
mild evidence for a triaxial shape, but the apparent triaxiality might 
be due to the presence of the Monocerus ring at low latitudes 
(|b| < 40°) and regions of the Southern Sagittarius stream. 

The root-mean-square deviation of the data around the max- 
imum likelihood model cr/tot typically ranges between 5% and 
20%. This indicates that the Milky Way stellar halo, or at least the 
component traced by the A-type stars in the SDSS DR8, is smooth 
and not dominated by unrelaxed substructure. 

This finding is discrepant with the conclusions of iBell et al.l 

( l2008h . These authors use main sequence turn-off stars selected 
from SDSS data release 5 (DR5) to model the steUar halo den- 
sity. They argue that no smooth model can describe the data and 
conclude that the stellar halo is dominated by substructure (with 
cr/tot ^ 0.33). We suggest that these contrasti ng results may be 
due to the different methods and tracers used. iBell et al.l ( |2008|) 
search for the lowest cr/tot fraction models. However, we find that 
the cr/tot values very little between different density models even 
if the likelihood is substantially increased. Furthermore, a lower 
cr/tot could also indicate that a model is 'fitting ' to any substruc - 
ture (e.g. triaxial models both in this work and in lBell et"al ]|2008h . 
iBell et al.l J200§) also use main sequence stars as tracers which are 
more numerous than A-type stars, but are much poorer distance in- 
dicators. The adopted absolute magnitude sc ale for main seque nce 
stars is dependent on metallicity and colour. iBell et al" as- 
sume a median absolute magnitude of Mr ^ 4.5 with a scatter of 
o'Mr = 0.9; an appropriate choice for a halo population with metal- 
licity [Fe/H] ~ —1.5. The scatter, even with an assumed metallic- 
ity, is greater than the spread in absolute magnitude for BS stars 
(cTA/g ~ 0.5 per colour bin). Uncertainties in the distances to stars 
could lead to compact substructures becoming more blurred ou t 
over a wider range of distance. It is possible that lBell et al.l | [200^ 
see a somewhat younger and/or metal-richer component of the stel- 
lar halo with the main sequence tracers. This component could be 
more unrelaxed than the one traced by A type stars. 

Our conclusion that the stellar halo is composed of a smooth 
underlying density, together with some additional substructures 
such as the Virgo Overdensity and the Sagittarius Stream, is very 
reassuring. If the stellar halo were merely a hotch-potch of tidal 
streams and unrelaxed substructures, then modelling and estima- 
tion of total mass and potential would be much more difficult. 
Many of the commonly-used tools of stellar dynamics - such as 
the steady-state Jeans equations - implicitly assume a well-mixed 
and smooth equilibrium. This raises the hope that a full understand- 
ing of the spatial and kinematic properties of stars in the smooth. 
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yet squashed and broken, stellar halo can yield the gravitational 
potential and dark matter profile of the Galaxy itself. 
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